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ABSTRACT 

Massive black hole binary coalescences are prime targets for space-based gravitational wave 
(GW) observatories such as LISA. GW measurements can localize the position of a coalescing 
binary on the sky to an ellipse with a major axis of a few tens of arcminutes to a few degrees, 
depending on source redshift, and a minor axis which is 2 — 4 times smaller. Neglecting weak 
gravitational lensing, the GWs would also determine the source's luminosity distance to better 
than percent accuracy for close sources, degrading to several percent for more distant sources. 
Weak lensing cannot, in fact, be neglected and is expected to limit the accuracy with which 
distances can be fixed to errors no less than a few percent. Assuming a well-measured cosmology, 
the source's redshift could be inferred with similar accuracy. GWs alone can thus pinpoint a 
binary to a three-dimensional "pixel" which can help guide searches for the hosts of these events. 
We examine the time evolution of this pixel, studying it at merger and at several intervals before 
merger. One day before merger, the major axis of the error ellipse is typically larger than its final 
value by a factor of ^ 1.5 — 6. The minor axis is larger by a factor of ^ 2 — 9, and, neglecting 
lensing, the error in the luminosity distance is larger by a factor of ^ 1.5 — 7. This large change 
over a short period of time is due to spin-induced precession, which is strongest in the final days 
before merger. The evolution is slower as we go back further in time. For z = 1, we find that 
GWs will localize a coalescing binary to within 10 deg^ as early as a month prior to merger 
and determine distance (and hence redshift) to several percent. 

Subject headings: black hole physics — galaxies: nuclei — gravitation — gravitational waves 



1. Introduction 

Among the most important sources of gravitational waves (GWs) in the low-frequency band of space- 
based detectors are the coalescences of massive black hole binaries (MBHBs) . Binaries containing black holes 
with masses in the range lO'' — 10^ Mq are predicted to form through the hierarchical growth of structure as 



dark matter halos (and the galaxies they host) repeatedly merge; see, for example, Sesana et al. (20071 and 



Micic et al. ( 2007[ ) for recent discussion. The Laser Interferometer Space Antenna {LISA) is being designed 



to have a sensitivity that would allow detailed measurement of the waves from these binaries. "Intrinsic" 
parameters — the masses and spins of the black holes which compose the binary — should be determined 
with very high accuracy, with relative errors typically ^ 10^"^ to 10^^, depending on system mass and 
redshift; see Lang & Hughes (2006, hereafter Paper I) for recent discussion. By measuring an ensemble of 
coalescences over a range of redshifts, MBHB GWs may serve as a kind of structure tracer, tracking the 
growth and spin evolution of black holes over cosmic time. 



"Extrinsic" system parameters, describing a binary's location and orientation relative to the detector, 
are also determined by measuring its GWs. In Paper I, we showed that a binary's position on the sky can be 
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localized at z = 1 to an ellipse with a major axis of a few tens of arcminutes and a minor axis a factor of 2 — 4 
smaller. At higher redshift (z = 3 — 5), these values degrade by a factor of a few, reaching a few degrees in 
the long direction and tens of arcminutes to a degree or two in the short one. We also found that, neglecting 
errors due to weak gravitational lensing, a source's luminosity distance typically can be determined to better 
1% at low redshift {z ~ 1), degrading to several percent at higher redshift (z = 3 — 5). 

The intrinsic ability of GWs to determine the distance to a coalescing binary is phenomenal. Coalescing 
MBHB systems constitute exquisitely well calibrated distance measures, with the calibration provided by 
general relativity. Unfortunately, this percent-level or better accuracy could only be achieved if we measured 
MBHB coalescences in an empty universe. In our universe, weak lensing will magnify or demagnify the GWs, 
and we will infer a luminosity distance smaller (for magnification) or larger (for demagnification) than the 
true value. This phenomenon affects all high-redshift standard candles. Its impact on Type la supernovae 



in particular has been discussed in detail ( ,Frieman 1997 Holz & Wald 1998 Holz 19981. It will not be 
possible to correct for this effect ( Dalai et al. 2003 1 since much of the lensing "noise" arises from structure 



on subarcminute scales that is not probed by shear maps (which map the distribution of matter on scales 
greater than 1' or so). Since we will not know the extent of the magnification when we measure MBHB 
waves, we must simply accept the fact that lensing introduces a dispersion of several percent in determining 
the distance to these GW events (see, e.g., Wang et al. [2002] to compute this dispersion as a function 
of redshift). When we quote distance measurement errors, we will typically quote only the intrinsic GW 
measurement error, neglecting lensing's impact. When the intrinsic GW distance error is < 5%, lensing will 
blur it to the several percent level. 

Note that a source's redshift z cannot be directly determined using only GWs. Gravitational wave 
measurements infer system parameters through their impact on certain dynamical timescales, such as orbital 
frequencies and the rate at which these frequencies evolve. Since these time scales all suffer cosmological 
redshift, z is degenerate with other parameters. For example, any mass parameter m is actually measured 
as (1 + z)m. However, if the binary's luminosity distance is determined, its redshift can then be inferred by 
assuming a cosmography. For most binaries, the redshift can be determined to several percent (with an error 
budget dominated by gravitational lensing. We thus expect that GW measurements will locate a binary 
to within a three-dimensional "GW pixel" which at z = 1 has a cross-sectional area of ~ 10~^ to 10~^ deg^ 
and a depth Az/z ~ several percent. 

It is anticipated that there will be great interest in searching the GW pixel for electromagnetic (e.g., 
optical. X-ray, radio) counterparts to MBHB GW events. Finding such a counterpart would be much easier 



if galactic activity were catalyzed in association with the coalescence (Kocsis et al. 20061. The nature of 



that activity is likely to depend rather strongly on the mass of the coalescing system ( Dotti et al. 2006 1 



For example, Armitage & Natarajan (2002) predict strong outflows and galactic activity prior to the final 



black hole merger as the smaller member of the binary drives gas onto the larger member, consistent with 
the high-mass (Mtot > lO'^M©) predictions of [Dotti et~al] (|2006|. |Milosavljevic fc Phinney| (|2005| describe 



an X-ray afterglow that would ignite after gas refills the volume that is swept clean by the coalescing binary; 
Dotti et al. predict this outcome for smaller systems (A/tot ^ several x 10^ Mq). Recent work by Bode & 



Phinney (20071 suggests that the final burst of radiation from a coalescing binary (which can convert ~ 10% 



of the system's mass to GWs very suddenly) may excite radial waves, and consequently electromagnetic 



^At redshifts z < 0.3, the error is actually dominated by peculiar velocity effects ( Kocsis et al.|[2006 l; however, the event 
rate is probably negligible at such low redshifts. As such, we will focus on gravitational lensing as the main source of systematic 
redshift error. 
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variability, in an accretion disk due to the quick change in the disk's Keplerian potential. Such a signature 
may be essentially mass independent. On the other hand, the coalescence may be electromagnetically quiet, 
in which case we face the potentially daunting task of searching the three-dimensional pixel for galaxies 
with morphology consistent with a (relatively) recent merger, or that have a central velocity dispersion a 
consistent with the inferred final black hole mass (assuming that the Mbh — o relation [Ferrarese & Merritt 
2000; Gebhardt et al. 2000] holds at the redshift of these sources, and so soon after merger). 

If the host galaxy or some other counterpart can be identified, we could then contemplate combining 
GW information with electromagnetic data. For instance, combining LISA mass measurements with the 



luminosity of the counterpart may allow us to directly measure the Eddington ratio L/L^dd (Kocsis et al. 



2006 ). Identifying a counterpart would also allow us to more accurately characterize the system. For example. 



much of the intrinsic luminosity distance error is due to correlations between distance and sky position. 
Finding an electromagnetic counterpart essentially determines a binary's position precisely, breaking those 
correlations. Previous studies have found that intrinsic distance error can be reduced by almost an order of 



magnitude if the position is known (Hughes 2002 Holz & Hughes 2005). (Lensing errors still dominate in 
such a case, so that the distance remains determined only at the few percent level.) A counterpart may also 
make it possible to simultaneously determine a source's luminosity distance and redshift. Such a "standard 
siren" (the GW analog of a standard candle) would very usefully complement other high-redshift standard 
candles ( |Holz fc Hughes 20051, such as Type la supernovae (iPhillips 1993 Riess et al. 1995 Wang et al. 
2003). A direct measurement of redshift will also break the mass-redshift degeneracy more accurately than 



can be done with just the luminosity distance and some assumed cosmological parameters. Breaking this 
degeneracy is critical when studying the growth of black holes with cosmic time ( |Hughes||2002| ). 

Many analyses ( |Cutler|[l998l |Hughes|[2002l |Vecchio|[2004l [Berti et al.|[2005l |Holz fc Hughesl[2005l ) have 
quantified how well LISA can determine MBHB parameters, including sky position and distance, using 
maximum likelihood Fisher matrix estimation. Our results from Paper I, given earlier, include the effects of 
"spin-induced precession" — precession of both the orbital plane and the individual spins of the black holes 
due to post-Newtonian spin-interaction effects. A significant result from that analysis is that spin- induced 
precession improves sky position accuracy by about half an order of magnitude in each direction versus 
previous analyses. This result can be partially understood as due to the breaking of a degeneracy between 
position and orientation angles: thanks to precession, the binary's orientation evolves with time and can be 



untangled from sky position. This effect was already known due to pioneering work by Vecchio (20041; by 



taking the analysis to higher order and considering a broader range of sources, we were able to show that 
this improvement held for essentially all astrophysically interesting MBHB sources. 

The purpose of this paper is to examine the localization of MBHB systems more thoroughly, in particular 
how the GW pixel evolves as the final merger is approached. Paper I only presented results for measurements 
that proceed all the way to merger. It will clearly be of some interest to monitor potential hosts for the binary 
event some time before the merger happens; if nothing else, telescopes will need prior warning to schedule 
observing campaigns. Understanding the rate at which localization evolves can also have an important 
impact on the design of the LISA mission, clarifying how often it will be necessary to downlink data about 
MBHB systems in order to effectively guide surveys. 

Our main goal is to understand for what range of masses and redshifts prior localization of a binary 
using GWs will be possible. A previous analysis by Kocsis et al. (2007b; hereafter K07) also examined this 
problem in great detail, but without including the impact of spin-induced precession. One of our goals is to 
see to what extent precession physics changes the conclusions of K07. We find that precession has a fairly 
small impact on the time evolution of the GW pixel except in the last few days before the final merger, at 
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which point its impact can be tremendous. Precession typically changes the area of the sky position error 
ellipse by a factor of ~ 3 — 10 (up to ~ 60 in extreme cases) in just the final day. This is in accord with the 
predictions of K07 (and even earlier predictions by N. Cornish 2005, unpublished). 

The structure of this paper is as follows. First, in § [2] we briefly review the basics of the MBHB 
gravitational waveform and the parameter estimation formalism that we use; this section is essentially a 
synopsis of relevant material from Paper I. Section |2.1| reviews the form of the GWs that we use in our 
analysis, while § |2.2| describes how those waves are measured by the LISA constellation. We describe the 
measurement formalism we use in § |2.3| In § |2.4[ we summarize our results from Paper I regarding the final 
localization accuracy that LISA can expect to achieve. 

We turn to a detailed discussion of the time evolution of the GW pixel in § |3] We begin by summarizing 
the key ideas behind the "harmonic mode decomposition" of K07 in § |3.1[ This technique cleverly allows 
calculation of the GW pixel and its time evolution with much less computational effort than our method 
(albeit without including the impact of spin precession). Unfortunately, we have discovered that some of the 
approximations used by K07 introduce a systematic underestimate of the final sky position error by a factor 
of 2 — 4 or more in angle; the approximations are much more reliable a week or more prior to the black holes' 
final merger. Modulo this underestimate, the K07 results agree well with a version of our code which does 
not include spin precession (particularly a week or more in advance of merger, when their underestimate is 
not severe). K07 thus serves as a useful point of comparison to establish the impact of precession on source 
localization. 



Section [3T2] is dedicated to our results, including comparison to K07 when appropriate. We find that all 
relevant parameter errors decrease slowly with time until the last day before merger, when they drop more 
dramatically. This sudden drop is not found in K07, nor is it present in a variant of our analysis that ignores 
spin precession. It clearly can be attributed to the impact of precession on the waveform. Before this last 
day, the major axis is 1.5 — 6 times, the minor axis ~ 2 — 9 times, and the intrinsic error in the luminosity 
distance Dj^ ~ 1.5 — 7 times bigger than at merger for most binaries (i.e., all except the highest masses). 
Going back to one week (one month) before merger, these numbers change to 2 — 9 (4—11) for the major 
axis, 3 — 14 (5 — 24) for the minor axis, and 3 — 14 (5 — 18) for the error in the luminosity distance. As a 
result, for z = 1, most binaries can be located within a few square degrees a week before merger and 10 deg^ 
a month before merger. The intrinsic distance errors are also small enough this early that Az/z remains 
dominated by gravitational lensing errors of several percent. Advanced localization of MBHB coalescences 
thus seems plausible for these binaries; the situation is less promising for sources at higher redshift. 

As a corollary to our study of the time evolution, we also examine the sky position dependence of errors 



(§ 3.3). The errors depend strongly on the polar angle with respect to the ecliptic, increasing in the ecliptic 
plane to as much as 35% over the median for the major axis, 85% over the median for the minor axis, and 
15% over the median for errors in the luminosity distance. The errors have a much weaker dependence on 
the azimuthal angle. When we convert to Galactic coordinates, we find that the best localization regions 
appear to lie fairly far out of the Galactic plane, offering hope that searches for counterparts will not be too 
badly impacted by foreground contamination. 

We conclude this paper in §[4] Besides summarizing our results, we discuss shortcomings of this analysis 
and future work which could help to better understand how well GWs can localize MBHB sources. 

Throughout the paper we set G = c = 1; a convenient conversion factor in this system is IO^Mq = 4.92 
s. When discussing results, we always quote masses as they would be measured in the rest frame of the 
source. These masses must be multiplied by a factor of 1 + z when used in any of the equations describing 
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the system's dynamics or its GWs (particularly the equations of §|2|. We convert between distance and 
redshift using a spatially flat cosmology with il\ = 0.75 (and hence ^Im ~ 0.25) and Hubble constant 
Ho = 75 km s"^ Mpc"^ 



2. Determining the sky position and distance accuracy of LISA 

In this section, we present necessary background. We begin by summarizing the inspiral GW signal and 
its measurement by LISA. We then describe the parameter estimation formalism that we use in our code. 
More details on these topics can be found in Paper I and references therein. Finally, we present a synopsis 
of the results from Paper I, summarizing the final sky position and distance accuracy that LISA can expect 
to achieve. Detailed discussion of how this accuracy evolves as a function of time is deferred to §|3] 



2.1. MBHB gravitational waves 



As is by now traditional, we break the GWs generated by a coalescing binary black hole system into 
three more or less distinct epochs: (1) a slowly evolving inspiral, in which the black holes gradually spiral 
toward each other as the orbit decays due to GW emission; (2) a loud merger, in which the black holes come 
together and form a single body; and (3) a ringdown, in which the merged remnant of the binary settles 
down to its final state. Our analysis focuses on the inspiral, the most long-lived epoch of coalescence and 
the epoch in which spin precession plays a major dynamical role. The inspiral's duration means that LISA's 
position and orientation changes considerably during the measurement; this, plus the impact of precession, 
is largely responsible for pinning down the sky position of coalescing sources. 



The inspiral GW signal is modeled using the post-Newtonian (PN) formalism ( Blanchet|2006 Blanchet 
et al. 1995 Will & Wiseman 1996). More specifically, we use the "restricted" PN approximation. The 



waveform can be written schematically as ( Cutler fc Flanagan||1994 ) 



h{t) = Re 



im*orb(*) 



(1) 



where x labels PN order, to is a harmonic index, and <i>orb(^) is the orbital phase. In the restricted PN case, 
we discard all amplitudes except the strongest Newtonian quadrupole, /ij. The phase, however, is computed 
to a desired PN order (here the second one, or 2PN). It has been recognized for some time that important 
additional information is carried by subleading harmonics (e.g., Hellings & Moore 2003). Recent work (Arun| 



et al.|[2007 Trias fc Sintes|[2008 1 has confirmed that sky position and distance are more sharply constrained 
when these terms are included in the wave model (although there still appears to be some disagreement 
among different groups as to the extent of the improvement). We plan to include these terms in future work, 
examining how higher harmonics and precession combine to localize MBHB systems. 

Spin-induced precession arises from geodetic and gravitomagnetic infiuences on the angular momentum 
vectors present in the binary. Averaged over an orbit, the leading-order spin precession equations are 
( [Apostolatos et al.||1994[ |Kidder|[T995 1 
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where mi and m2 are the masses of holes 1 and 2, M = mi + m2 is the total mass, /i = mim^/M is the 
reduced mass, L is the normal to the orbital plane. Si and S2 arc the spins of the holes, and r is their 
orbital separation in harmonic coordinates. Since total angular momentum is conserved on short timescales 
(ignoring the inspiral), the orbital angular momentum precesses to compensate. These precessions drive 
several modifications to the orbital phase. One correction is due to spin-orbit and spin-spin terms in the 
post-Newtonian expression for orbital frequency: Precession makes these terms time dependent, adding a 
new modulation to the GW phase. Another correction arises from the changing orientation of the orbital 
plane ( Apostolatos et alT]|1994 l. 



The waveform is described as two polarizations propagating in the — n direction, where n gives the 
direction to the binary on the sky: 

2A/f5/3C7rf^2/3 

h+{t) = ^J' [l + {L-hr]cosmt) + Sr>m], (4) 

hx{t) = - n ^ (L.n)sin[<i>(t)+ W)]. (5) 

where / is the frequency of the waves, Dl is the luminosity distance, ^{t) = J* 2TTf{t')dt' is the wave phase 
(twice the orbital phase in the quadrupole approximation), and Sp^{t) is the precessional correction to the 
orbital phase due to the changing orientation of the orbital plane. The "chirp mass" A4, which controls the 
frequency evolution, or chirp, of the signal, is given by = ir'/^Np/^ . Note that the relative weight of the 
polarizations is determined by the binary's inclination to the line of sight t, where cos t = L • n. 



2.2. MBHB gravitational waves as measured by LISA 

The baseline LISA mission consists of three spacecraft arranged in an equilateral triangle with arms of 
length L — bx 10® km. The center of this constellation orbits the Sun 20° behind the earth, with the triangle 
oriented at 60° to the ecliptic. The individual spacecraft thus orbit in different planes, causing the triangle 
to "roll" about itself in its orbit. To model LISA's response to a signal, we use the long wavelength (A ^ L) 



approximation introduced by Cutler ( 1998 1. The three arms of the triangle act as a pair of two-arm detectors; 
the signals from the three arms are combined in such a way that the noise spectrum in each "synthesized" 
detector is independent of the other. Each detector measures both GW polarizations, weighted by antenna 
response functions F^[t) and F^{t), where i G {1,11} labels the synthesized detector. These functions 
depend on the source's sky position as measured in the detector fTaine. Because of the constellation's rolling 
orbital motion, these functions are time dependent. The motion around the Sun also causes a Doppler shift 
in the signal's frequency. This time dependence allows the sky position of the source to be measured. 

An additional time dependence is produced by the precession of the binary's orbital plane. We saw above 
that the polarizations strongly depend on the binary's inclination to the line of sight l (via cos t = L • n) . 
The antenna pattern functions F^'^ (t) also depend on the component of L that is perpendicular to the line 
of sight. Since precession changes the direction of the orbital plane, the measured polarizations would vary 
in time even if the detector were static. The timescale of these precessions is much shorter than the LISA 
orbital timescale of 1 yr, so they provide a great deal of additional information about sky position. 



We perform the parameter estimation analysis in the frequency domain, so we Fourier transform the 
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signal using the stationary phase approximation. The final measured waveform is given by 

where ^poi,i(i) is the "polarization amplitude," <^poi,i(^) is the "polarization phase," and ipoit) is the Doppler 
phase. The first two encode modulations due to the relative polarization weightings and antenna pattern 
functions; the latter is due to the Doppler shift in the GW frequency. Full expressions for these functions 
can be found in Paper I. To second post-Newtonian order, the phase ^'(/) is given by dPoisson fc Will|1995| ) 



20 /743 11 ^ , ,,,,2/3 



,„/3058673 5429 617 , \, ,^.^4/^ 



(7) 



where -q = /x/M, is the "merger time" (the time at which / formally diverges in the post-Newtonian 
framework), <i>c is the phase at tc, and (3 and a are "spin-orbit" and "spin-spin" parameters, respectively, 
encoding the manner in which interactions between the spins and the orbits impact the phase. See Paper I 
for precise definitions. For our purposes, the most important property of (3 and a is that they oscillate due 
to precession, modulating the phase ^'(/). Note that waveform Q depends on 15 parameters: 2 masses, 6 
initial spin components, 1 distance, 2 sky angles, 2 initial orbital plane angles, and <i>c- 

The end of inspiral/beginning of merger is a somewhat ad hoc and fuzzy boundary. Indeed, recent 
numerical computations have shown that the GWs produced by a binary that coalesces into a single body do 
not show any particular special feature as the black holes come together, instead smoothly chirping through 
this transition ( Buonanno et aL]|2007 Pan et ar]|2008 1. Since we are not including numerical merger waves 



in our analysis, we require some point to terminate our post-Newtonian expansion. Most studies show that 
the inspiral comes to an end when the separation of the bodies in harmonic coordinates is roughly r ~ 6M; 
at this point, the system's GW frequency is approximately given by 

2 0.02 

/merge - ;^^^Kopler(f' = 6M) ~ — - , (8) 
ZTT IVl 

where f^Kepior — {M/r^Y^'^ is the Keplerian orbital angular frequency. (The factor 1/2. converts from 
angular frequency to frequency; the additional factor of 2 accounts for the quadrupolar nature of gravitational 
waves.) We use equation (pi) throughout our analysis to terminate the inspiral. 



2.3. Summary of parameter measurement formalism 



To estimate how well LISA will be able to measure MBHB parameters, we use the maximum likelihood 
formahsm developed by Finn ( 1992[ ). The GW signal hi{t) described above is corrupted by noise to produce 
a measured signal Siif) = hi(t) + ni{t). We assume the noise ni(t) is zero mean, wide-sense stationary, 
Gaussian, and uncorrelated between the two detectors (i G {I, II}). It can be characterized by its (one 
sided) power spectral density (PSD) S'„(/), which is described in detail in § IIIB of Paper I. The PSD 



includes instrumental noise intrinsic to the detector, based on the model described by Larson et al. (20001 



and implemented on the World Wide WelJ^ with important numerical factors described in § IIC of Berti 



2 See 



http;/ /www. srl.caltech.edu/ shane/sensitivity/ 



et al. 



et al 



( 2005 1 . It also includes astrophysical noise arising from a confused background of Galactic ( Nelemans 



20011 and extragalactic ( Farmer fc Phinney||2003 1 white dwarf binaries. With the noise PSD, we can 



define the inner product between two signals a and b as 



{a\b) =4 Re 



df 



a*imf) 

Snif) 



df 



a*imf)+&if)b*if) 
Snif) 



(9) 



We write a GW as h{9), where the components of the vector 9 represent the 15 parameters describing 
MBHB waves. We assume that a detection has occurred — that a sign al with particular parame ters is 
present in the data — and that maximum likelihood (ML) estimates 6 (Cutler & Flanagan 19941 of these 
parameters have been found. The signal-to-noise ratio (S/N) is then given by 



1/2 



(10) 



Our ultimate goal is to understand the errors in the ML parameters. To do so, we examine the probability 



that the GW parameters are given by 6, expanded around the ML estimate 6 (Cutler & Flanagan 
Poisson fc W"ml[l995l ): 



1994 



p(0|s)cxexp(--r„65r^ 



where (56'" ^ 9" - 9" and 



dh 



W J ' 



(11) 



(12) 



evaluated a,t 6 — 6, is known as the Fisher matrix. This expression holds for large values of the S/N and is 
known as the Gaussian approximation. For the two-detector case, we can exploit the independence of the 
noises to write the total Fisher matrix F*"^* as the sum of the individual Fisher matrices. If we define the 
covariance matrix S"'' = (Fj^iJ:)'''', we can then estimate the measurement error in some parameter 0", 



A61" = ^/{{6¥y 



(13) 



as well as the correlation between two parameters. 



A9°-A9'' y/j^aaj^bb 



(14) 



We are particularly interested in the sky position, which we parameterize by the coordinates {p,N = 
cos^AT, 0Ar)j^ We want to convert from errors in these two parameters to an error ellipse on the sky. 
To do so, we first perform a change of coordinates from p,N to 9^. For small deviations from the ML 
estimate, we have 69n — id9N /dfiN)6fLN = — (5/ijv/sin07v- Next, we recognize that due to the geometric 
properties of the sphere, the same 64)n corresponds to a different "proper" angle depending on the value of 
9^: 6(f>^ — sin0jv^0Af- With these modifications, equation (111 becomes 



p{e\s)^e^pf-^-Tl,,,69;'69: 



(15) 



^The bar over the angles conforms to the notation in Paper I for quantities measured in the solar system barycenter frame, 
as opposed to those measured in the time- varying detector frame. 
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Here (50° = {69 n, S(f>^) denotes the proper errors accounting for the metric of the sphere, and T^if^t represents 
the equivalent Fisher matrix with all conversion factors absorbed inside: 



4'N<t>N 



= CSC'' 9 nT 



- rP 



= (-sin^Ar)(cSC^Ar)rp^ 



(16) 
(17) 
(18) 



and so on for the rest of the elements. The inverse of this matrix is the proper covariance matrix, Sp . 
Consider now just the 2x2 subspace of the covariance matrix containing the sky position variables. Let the 
eigenvalues of this subspace be A±. If we define the error ellipse such that the probability that the source 
lies outside of it is e^^ (corresponding to a « 63% confidence interval), then the major and minor axes are 
given by 2a — 2(2A_|_)^/^ and 26 = 2(2A_)^/^, respectively. Expressed in terms of the original covariance 
matrix, these are 



csc^ + sin^ ± d {csc^ e^Et^^f^^ - sin^ 5)^5]0iv0iv)2 _^ 4(5]MJv0iv)2 



1/2 



(19) 



Many previous analyses have reported the ellipse's area Afijv = 7ra& or (Afi^v)^''^, the side of a square of 
equivalent the sky position error ( |Cutler|1998l|Vecchio|2004l|Berti et al.|2005| [HoTz fc Hughes|2"005| . 

Information about the ellipse's shape, crucial input to coordinating GW observations with telescopes, is not 
included in such a measure. By examining both 2a and 26, this information is restored. 

In an empty universe, the third dimension of our pixel would be determined from the luminosity dis- 
tance error, AD^/Dl = (5]in_Di,,inDi,-|i/2^ already discussed, this quantity is often so small that weak 
gravitational lensing is expected to dominate the distance error budget. We thus expect 



Dl 



(20) 



where the lensing dispersion crions(-2) is described in, for example, Wang et al. (20021. If cosmological 
parameter errors can be neglected, then one typically finds that the redshift error is about equal to the 
distance error (|Hughes||2002|): Az/z « ADl/Dl, independent of redshift. 



2.4. Final position and distance accuracy 

We conclude this section by reviewing the final accuracy with which MBHB sky position and distance can 
be determined with GWs. In Paper I, we described the code which implements the measurement formalism 
described in §§ |2.1| - 12.3| above. This formalism is coupled to a Monte Carlo engine: We specify the rest 
frame masses and redshift and then randomly choose the sky position, initial angular momentum and spin 
directions, and spin magnitudes for each binary. We also uniformly distribute the time parameter tc of 
each binary over the assumed duration of the LISA mission (which we take to be 3 yr). Some binaries are 
therefore already "on" when the mission begins, and as such may radiate in band for less time than other 
binaries which enter the band later. The random distribution of merger times also means that the relative 
azimuth between a binary's sky position and LISA^s orbital position at merger, Scj) — (j)N — 4>LisA{ic), is itself 
randomly distributed, even if we choose to focus on a particular ^^r. We discuss the implications of this in 



more detail in § 3.3 Computing the errors for a large number of binaries (10* for each mass and redshift 
set, unless otherwise stated) allows us to make statements about the distribution of parameter errors for 
different mass and redshift combinations. 
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Table [T] summarizes the final accuracy of MBHB sky position and distance measurements. Here we 
show the median measurement accuracy from many Monte Carlo surveys. Several trends are particularly 
noteworthy. At all redshifts, the accuracy is worst for the largest masses. This is because the most massive 
systems are in band for the least amount of time: The frequency at which inspiral ends is inversely pro- 
portional to mass (see eq. [s]), and more massive systems move more rapidly from low frequency to high 
frequency. The short time these systems spend in band means that LISA measures a relatively small number 
of modulations (whether induced by the constellation's orbit or by spin precession). Second, note that the 
results for mi — m2 tend to be less accurate than results with similar total mass but for which mi > 7712. 
The cause of this phenomenon lies in precession equations ([2| and ^ : When mi 7^ m2 , the two spins precess 
at different rates, imposing richer modulations on the measured GWs. Since mi ~ m2 is a rather artificial 
limit, we expect that the more accurate results for nonunity mass ratio will be the rule. 

Independent of these trends, an important result is that MBHB systems are pinned down on the sky fairly 
accurately at z = 1. Modulo the higher mass binaries, the median major axis of the sky position error ellipse 
is typically about 15' — 45', and the median minor axis is about 5' — 20', with a total ellipse area considerably 
smaller than 1 deg^ (ranging from about 0.02 to 0.2 deg^). Sources at z = 1 are located accurately enough 
that one can comfortably contemplate searching the GW error ellipse for MBHB counterparts with future 
survey instruments, such as the Large Synoptic Survey Telescope (LSST; Tyson et al. 2002). 

At higher redshift, positional accuracy degrades. This is due to the weakening of the signal with distance 
and to the redshifting of the waves' spectrum, so that the signal tends to spend less time in band. At z = 3, 
the major axis of the error ellipse is ~ 1° — 4° across, and the minor axis is ~ 40' — 110'. The total area 
of this ellipse is ~ 0.5 — 5 deg'^. At 2 = 5, this degrades further to ~ 3° — 5° for the major axis, ~ 1° — 3° 
for the minor axis, and ^2 — 10 deg^ for the total area. These degraded accuracies are still sufficiently well 
determined that telescopic searches for MBHB counterparts have a good chance to be fruitful (although not 
nearly as simple as they would be at z ^ 1). 

In all cases, the GW distance determination is extremely accurate: For all but the highest masses, 
SDl/Dl ^ 0.7% at z = 1, < 4% at z = 3, and < 5% at z = 5. Distance is determined so precisely that 
these errors are in fact irrelevant — weak gravitational lensing will dominate the distance error budget for 
all but the most massive MBHB events. 

Although the median values reported in Table [l] indicate the typical sky position accuracies we expect, 
it should be emphasized that they are taken from broad distributions. Figure [T] presents the distributions 
we computed for binaries at z = 1 with masses mi = 10^ M0 and m2 — (10^,3 x 1O^,1O^)M0. Note that 
the major axis distribution {left) is rather flat when compared to the minor axis distribution (right). It lacks 
a single well-deflned peak; in fact, it is actually bimodal for mi/m2 > 1, with one peak near 10' — 20' and 
another closer to 1° — 2°. We find that this behavior holds over all mass and redshift cases of interest, with 
only slight variations. At smaller masses, the distribution is broader than the cases pictured, without strong 
bimodality; for larger masses, the distribution is somewhat narrower and tends to develop two very distinct 
peaks. Higher mass ratios tend to accentuate the peaks. These results hold for higher redshift, except that 
transitions between various behaviors occur at smaller total mass [in keeping with the fact that it is not 
mass but (1 + z) times the mass that determines dynamical behavior]. 

The minor axis distribution exhibits a rather long tail to very small values, especially when the mass 
ratio is large. For mi/m2 — 10, the distribution peaks near a minor axis ^ 10', but extends from roughly 
1" to about 100'. As the mass ratio approaches 1, the peak moves to slightly larger values (slightly more 
than 10' for mi/m2 = 3; roughly 30' for mi/m2 = 1), and the tail becomes less populated (although the 
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Table 1: Sky position and distance accuracy at merger 



z 


TOi (Mq) 


mi (Mq) 


2a (arcmin) 


26 (arcmin) 


ArjAr(dcg2) 




1 


10^ 


10^ 


27.3 


13.3 


0.0729 


0.00398 




3 X 10^ 


10^ 


16.9 


7.33 


0.0233 


0.00240 






3 X 10^ 


23.3 


11.8 


0.0556 


0.00357 




10^ 


10'' 


27.2 


6.62 


0.0235 


0.00320 






3 X 10^ 


31.3 


13.2 


0.0705 


0.00393 






10^ 


40.2 


21.9 


0.176 


0.00560 




3 X 10^ 


3 X 10^ 


34.1 


9.20 


0.0445 


0.00376 






10*^ 


32.3 


14.7 


0.0839 


0.00419 






3 X 10^ 


43.3 


22.3 


0.193 


0.00689 




10^ 


10^ 


37.6 


12.2 


0.0670 


0.00457 






3 X 10^ 


42.1 


19.0 


0.142 


0.00610 






10^ 


81.3 


38.6 


0.680 


0.0250 


3 


10^ 


10^ 


81.0 


40.8 


0.665 


0.0123 




3 X 10^ 


10^ 


92.5 


39.5 


0.656 


0.0126 






3 X 10^ 


142 


75.7 


2.15 


0.0201 




10^ 


10^ 


141 


36.6 


0.739 


0.0155 






3 X 10^ 


129 


56.7 


1.25 


0.0161 






10^ 


158 


84.3 


2.64 


0.0237 




3 X 10^ 


3 X 10^ 


132 


40.3 


0.751 


0.0153 






10^ 


142 


64.6 


1.65 


0.0193 






3 X 10^ 


224 


111 


5.08 


0.0422 




10' 


10^ 


206 


78.5 


2.74 


0.0293 






3 X 10® 


297 


152 


9.40 


0.0805 






10^ 


2000 


583 


256 


2.41 


5 


10^ 


10^ 


169 


85.7 


2.93 


0.0260 




3 X 10^ 


iqs 


217 


95.8 


3.73 


0.0284 






3 X 10^ 


295 


161 


9.29 


0.0409 




lO'^ 


10^ 


248 


66.8 


2.35 


0.0273 






3 X 10^ 


233 


101 


3.96 


0.0294 






10® 


315 


162 


10.2 


0.0501 




3 X lO^ 


3 X 10^ 


265 


86.4 


3.27 


0.0318 






10*^ 


304 


139 


7.52 


0.0436 






3 X 10*^ 


538 


260 


29.5 


0.140 




10^ 


10® 


577 


290 


31.9 


0.124 






3 X 10® 


1720 


621 


234 


1.24 






10^ 


180000 


29600 


1.15 X 10® 


377 



Note. — Median errors in sky position and distance for binaries at various masses and redshifts z = 1, 3, and 5: '2a is the 
major axis of the sky position error ellipse, 26 is the minor axis, and Af2/\? is the elhpse's area. (Note that since all of these 
quantities are medians, it is not true that AQjv = 7ra6.) At each mass and redshift point, 10'' binaries are drawn with randomly 
distributed sky positions, initial orbit and spin orientations, spin magnitudes, and tc- 
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Fig. 1. — Distribution of the major axis 2a (Ze/<) and minor axis 26 {right) of the sky position error eUipse 
for lO'* binaries with mi = 10^ Mq and m2 = 10^ Mq {solid line), 3 x 10^ Mq {dashed line), and 10^ Mq 
{dash-dotted line) at z = 1. Note the bimodal character of the major axis as well as the long tail down to 
small minor axis, both of which are particularly prominent for larger mass ratios. 



distributions span roughly the same extent as when r7ii/m2 — 10). We find that this tail exists for all 
interesting mass and redshift combinations, with the same strong dependence on mass ratio as shown in 
Figure [T] 

Finally, it is worth emphasizing that spin-induced precession has a significant effect on the accuracy 
which we report here. In Paper I, we compared these results to those obtained with a code which does not 
include precession physics. Position and distance accuracy are determined only by the detector's motion in 
this case. We found that precession effects reduce the major axis of the sky position error ellipse by a factor 
of ^ 2 — 7 and the minor axis by a factor of ~ 2 — 10. The distance error is likewise improved by a factor of 
~ 2 — 7. Factors of a few or more improvement are also seen when additional wave harmonics are included 
( Arun et "ar]|2007 Trias fc Sintes|2068 1 . We are eager to see whether the improvements from precession and 
from higher harmonics can be combined. 



3. Time and location dependence of the GW pixel 

We turn now to a detailed discussion of how well GWs localize an MBHB system as a function of 
time before final merger and as a function of sky location. We begin by discussing the analysis of K07, 
which presents a clever algorithm for estimating extrinsic parameter errors as a function of time until merger 
(although at present it does not include spin precession). We demonstrate that their analysis unfortunately 
underestimates final position errors by roughly a factor of ~ 2 — 4 or more (in angle) due to neglect of 
certain parameter correlations; the underestimate is much less severe a week or more prior to merger. We 
then present our own results for the time dependence of the GW pixel, using K07 for comparison where 
appropriate. Finally, we conclude with a brief study of the pixel's dependence on sky location. 
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3.1. Summary of K07 

K07 have devised a new method, the harmonic mode decomposition (HMD), to solve for the extrinsic 
parameter errors as a function of time to merger. In the HMD, modulations caused by LISA^s motion 
are decoupled from the much faster inspiral timescales and are then expanded in a Fourier series. The 
resulting expression for the measured signal features a time-dependent piece with no dependence on the 
extrinsic parameters and a time-independent piece with all the parameter dependence. As a result, when 
Monte Carlo simulations are done across parameter space, the time-dependent integrals do not need to be 
recomputed for each sample of the distribution. This makes it possible to quickly survey the estimated 
parameter errors across a wide range of parameter space. 

As already emphasized, the waveform model used by K07 does not (yet) include the impact of spin 
precession. As such, we intend to use their results as a baseline against which the impact of spin precession 
can be compared. Before doing so, we first checked to make sure that their results were in agreement with 
a variant of our code which does not include spin precession fHugh es|[2002| . To our surprise, we found that 



the final position accuracy predicted by K07 was typically a factor ~ 2 (in angle) more accurate than our 
code predicted. 

After detailed study of the HMD algorithm and comparison with our (precession free) code, we believe 
we understand the primary source of this disagreement. K07 define a set of "slow" parameters ^siow, which 
correspond (with some remappings) to our extrinsic parameters: IhDl, cosOl, cosOn, 4>Li and (/)jv. Here 
i^Li 't'L) defines the orientation of the orbital angular momentum L, which is constant when spin precession 
is ignored. K07 also define a set of "fast" parameters ©fast, which, modulo the exclusion of spin, map to 
our intrinsic parameters: tc, ^c, InAl, and Inry. In their formulation of the HMD, K07 approximate the 
cross-correlation between intrinsic and extrinsic parameters to be zero. Although the correlations between 
intrinsic and extrinsic parameters tend to be small, they are not zero. We find that they typically range in 
magnitude from about 0.1 to 0.4, sometimes reaching ~ 0.8. Neglecting these correlations altogether leads 
to a systematic underestimate in extrinsic parameter errors. 

An example of this is shown in Tables |2] and [3] To produce the data shown in Table [2] we compute the 
Fisher matrix T'^°^ and then invert for the covariance matrix, E"^ = {T^^tf^. Table^then presents a slightly 
massaged representation of this matrix: Diagonal elements are the 1 a errors (S"")^/^, and off-diagonal 
elements are the corr elation coefficients c'^^ = S'^''(E''''E^'')-i/^ Take particular note of the magnitude of 
the correlations between intrinsic and extrinsic parameters (the upper right-hand portion of Table |2]). Many 
entries have values ~ 0.2 — 0.3, and two have values ~ 0.7 — 0.8. 

To see what effect neglecting the intrinsic-extrinsic correlations has, we repeat this exercise, with a 
slight modification: We compute F*"^* as before, but we now set to zero entries corresponding to mixed 
intrinsic/extrinsic parameters. For example, we set by hand Finu^, ii\M — 0. We then invert this matrix 
to obtain S"''. The result is shown in Table |3] Note that mean parameter error (diagonal entries) is often 
significantly smaller than errors when these correlations are not ignored. The impact of correlations between 
intrinsic and extrinsic parameters is clearly not negligible. 

Table [4]gives further examples illustrating the impact of neglecting these correlations on our estimates of 
LISA'S localization accuracy. We show 10 points drawn from a 10^ binary Monte Carlo run; all use the same 
masses and redshifts (toi = 3 x 10® Mq, m2 = 10® Mq, and z = 1) but have different (randomly distributed) 
sky positions, orientations, spins, and tc- For these parameters, we find that neglecting intrinsic-extrinsic 
correlations causes one to underestimate the major axis of the position ellipse by a (median) factor of ~ 2 
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Table 2: Full errors and correlations 





In Dl 


cos 9l 


cos On 




0Ar 


tc 




InM 


In rj 




a 




0.233 


-0.984 


0.878 


0.509 


0.213 


0.0801 


0.246 


0.227 


-0.186 


0.205 


-0.106 


cos 9l 




0.467 


-0.861 


-0.350 


-0.071 


-0.040 


-0.098 


-0.178 


0.138 


-0.156 


0.0622 


cos On 






0.0006 


0.465 


0.203 


0.0709 


0.231 


0.201 


-0.166 


0.181 


-0.095 


h 








0.687 


0.782 


0.232 


0.827 


0.337 


-0.317 


0.328 


-0.259 


4>N 










0.0017 


0.193 


0.691 


0.252 


-0.244 


0.250 


-0.210 


tc 












63.1 


0.705 


0.923 


-0.955 


0.942 


-0.993 


$c 














4.00 


0.742 


-0.747 


0.748 


-0.726 


InM 
















0.0010 


-0.995 


0.998 


-0.956 


In 77 


















0.303 


-0.999 


0.981 


/9 




















1.11 


-0.971 


(T 






















0.722 



Note. — Example of errors (diagonal elements) and correlations (off-diagonal elements) for a binary with rrii = 3 X 10® Mq 
and 7712 = 10® Mq at 2 = 1. The errors in i^^, tpff, and <I'c are measured in radians; the error in tc is measured in seconds. 
This example was taken from the same Monte Carlo distribution used to make Table |4] in this particular case, the randomly 
distributed parameters have the values cosdj^ = —0.628, cosOpf = 0.850, 4>L = 3.50 rad, 4>N = 0.514 rad, tc = 6.90 X 10^ s, 
/3 = 1.48, and a = 0.107. Entries containing ellipses can be found by symmetry. 



Table 3: Errors and correlations neglecting correlations between intrinsic and extrinsic parameters 





\nDL 


cosOl cosOn 4>l 


0Af 


tc 




InM 


In rj 




a 


In Dl 


0.142 


-0.999 0.721 0.999 


0.0602 




















cos Ol 




0.294 -0.721 -0.999 


-0.0611 




















cos On 




0.0004 0.719 


0.0576 
























0.120 


0.0664 




















4>N 






0.0010 




















tc 








61.3 


0.990 


0.929 


-0.959 


0.946 


-0.993 


<^c 










2.10 


0.969 


-0.988 


0.981 


-0.999 


\nM 












0.0010 


-0.995 


0.998 


-0.960 


In r/ 














0.288 


-0.999 


0.983 


P 
















1.05 


-0.974 


a 


















0.697 


Note. — 


Exampl 


e of errors (diagonal elements) and correlations (off-diagonal elements) for a 


binary with 


mi = 3 X 10® Mq 


and m2 = 


10® Mq 


at 2: = 1 if correlations between intrinsic and extrinsic paramet 


srs are neg 


lected. The 


errors in 


4'L, 4'N, 



and <I>c are measured in radians; the error in tc is measured in seconds. This example was taken from the same Monte Carlo 
distribution used to make Table]!] in this particular case, the randomly distributed parameters have the values cos di^ = —0.628, 
cosSpf = 0.850, <j>L = 3.50 rad, ^jv = 0.514 rad, tc = 6.90 X 10^ s, f3 = 1.48, and a = 0.107. Entries containing ellipses can be 
found by symmetry. 



- 15 - 



Table 4: Example sky position error measures from Monte Carlo sample, comparing full Fisher matrix 
technique with the K07 approximation 



2a (arcmin) 


2b (arcmin) 


AO 


[acg ) 


Full 


K07 


Full 


K07 


Full 


K07 


201 


63.5 


59.0 


15.1 


2.59 


0.210 


165 


120 


108 


88.6 


3.90 


2.32 


117 


61.6 


81.0 


14.3 


2.07 


0.193 


197 


69.7 


46.7 


13.4 


2.01 


0.204 


10.9 


7.23 


8.21 


5.03 


0.0196 


0.00793 


197 


51.2 


55.2 


13.3 


2.37 


0.149 


46.7 


26.8 


36.5 


9.19 


0.372 


0.0538 


18.1 


12.7 


15.0 


6.37 


0.0595 


0.0177 


155 


92.4 


88.5 


16.2 


2.98 


0.326 


146 


143 


139 


10.9 


4.43 


0.342 



Note. — Ten Monte Carlo points comparing sky position accuracy for binaries with mi = 3 X 10^ Mq and m2 = 10^ A/q at 

2=1. 

and the minor axis by a factor of ~ 3 — 4; the area is underestimated by a factor of ~ 6 — 7. 

Our conclusion is that the HMD technique developed by K07 is overly optimistic by a factor of ^ 2 — 4 
or more (in angle) regarding the final accuracy with which GWs can locate an MBHB event on the sky. As 
a prelude to the time evolution study we present in § |3.2[ we also examined how this underestimate evolves 
as merger is approached. To our relief, it appears that this underestimate is much smaller prior to merger: 
For the handful of cases we examined, the factor of 2 — 4 underestimate in angle falls to a mere 10% — 25% 
offset one week prior to merger. We find that the offset plateaus at this level, remaining at a few tens of 
percent up to 28 days before merger. 

Accounting for this systematic final underestimate, we thus find K07's results to be a good baseline 
against which to compare our results. This comparison makes it possible to assess the extent to which spin 
precession improves our ability to locate massive black hole binaries prior to the final merger. 

3.2. Results I: Time evolution of localization accuracy 

We finally come to the main results of this paper, the time evolution of our ability to localize MBHB 
systems using GWs when spin precession is included. The results summarized in § |2.4| describe the size of the 
GW pixel (sky position error ellipse and luminosity distance) at the end of inspiral. We do not at this time 
incorporate any information regarding the merger and ringdown phases. The end-of-inspiraQ localization 
accuracies are good enough that searching the GW pixel for counterparts to MBHB coalescences is likely to 
be fruitful. However, given how little is understood about electromagnetic counterparts to these events, it is 



^Throughout this paper, "final" accuracy refers to the end of inspiral. 
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unclear if waiting until these final moments is the best strategy for such a search. It will surely be desirable 
to also monitor the best-guess location some days or weeks in advance for electromagnetic precursors to the 
final merger. The rate at which the GW pixel evolves as we approach the merger will have strong implications 
for determining the rate at which LISA data is sent to the ground. 

To this end, we now examine the time dependence of the LISA pixel. We have modified the code 
from Paper I to stop the calculation at a specified time before the fiducial merge frequency, equation (|8|. 
We begin the evolution of each binary at the moment it enters the LISA band (which we take to occur at 
/low = 3 X 10~^ Hz). Because we randomly distribute tc over our (assumed) 3 yr LISA mission, some sources 
are already in band at the mission's start; consequently, these sources begin at / > /low The binary's 
evolution is then followed until it reaches a GW frequency /stop = /(^(/merge) — N), where N is the number 
of days before merger that we want to stop the signal. (Choosing iV = duplicates the analysis of Paper I.) 

Figure [2] shows the error ellipse evolution for nine examples taken from a sample of 10^ computed with 
mi = 10*5 Mq, m2 = 3 X lO^M© and z = 1. We show results for = 0, 1, 2, 4, 7, 14, 21, and 28. For each 
binary, the major axis is plotted on the x-axis, while the minor axis is plotted on the j/-axis; we do not show 
how each ellipse would be oriented on the sky. The results shown in Figures |2]3 and|2^ were selected by hand 
from the distribution as examples of contrasting behavior. The binary in Figure [2]a shows a dramatic change 
in the error ellipse with time, especially in the last day before merger. In that day, the binary is localized 
to an ellipse with 2a — 6.67' and 2b = 6.25', an area ^ 60 times smaller than at = 1. By contrast, the 
binary in Figure |2^ shows almost no change in the error ellipse over the entire four weeks prior to merger. 

These are clearly extreme cases. Other extremes exist, including binaries with a minor axis orders of 
magnitude smaller than the major axis (see the tail in Fig. [l| right), binaries where the evolution of one or 
both axes is not strictly monotonic, and binaries which have very large ellipses (essentially filling the sky) 
for large N. (Such cases correspond to binaries which are already well into the LISA band when the mission 
starts; for large N, there is little baseline for the various modulations to encode their position.) To get a 
sense of more typical behavior, we selected the binaries in Figures [2]c - [2]! randomly from the distribution. 
There does not appear to be any "typical" evolution; each binary exhibits some unique features. Most 
binaries, however, seem to share with the binary in Figure |2]a the property that the final day before merger 
gives much more information on the position than any day before it (albeit to a lesser degree). We will see 
below that this feature holds for the medians of almost all mass and redshift cases. It is worth noting that 
although K07 agree with us on most of the other qualitative features of the time dependence, they do not 
see the dramatic change in the final day of inspiral. As we will discuss in more detail below, this dramatic 
improvement toward the end of inspiral is due to spin precession physics. This is in good agreement with 
the expectations of N. Cornish (2005, unpublished) and K07 that spin precession would most dramatically 
impact the last week or so of inspiral. 

While the evolution of parameter errors for individual binaries is interesting, of more relevance is the 
evolution of the errors' distribution. The left panel of Figure [3] shows the time dependence of the distribution 
of the major axis 2a for our model system of nii = IO^AIq, m2 = 3 x lO^M©, and z = 1. We can clearly 
see the evolution to smaller major axis as the binary nears merger. Four weeks before merger, the median 
major axis is 4.8 times larger than at merger; this number shrinks to 3.9 two weeks before merger, 3.2 one 
week before, and 2.5 two days before. As expected from the individual binaries, the most dramatic change in 
the distribution occurs during the last day before merger. Not only is the median substantially reduced (by 
a factor of 2.2), but the shape sharply changes. For > 0, the distribution is distinctly peaked, becoming 
gradually flatter as N gets smaller. Over the last day of inspiral, the distribution evolves into the almost 
entirely flat, slightly bimodal shape first seen in Figure [T] We find that this same behavior holds for all masses 
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Fig. 2. — Evolution of the sky position error ellipse for nine individual binaries selected from a set of 10"*. All 
have TOi = 10^ Mq, m2 = 3 x 10^ Mq, and z = 1. The ellipses are oriented so their major axes are parallel to 
the a;-axis and their minor axes are parallel to the y-axis; the axes are labeled in arcminutes. From outside 
in, the ellipses are evaluated at 28, 21, 14, 7, 4, 2, 1, and days before merger. 



and redshift cases of interest: A sharply peaked distribution at > evolves into the flatter, sometimes 
bimodal final distributions described in § |2.4[ As the total mass increases, however, the final distributions 
become so narrow that the shape change is no longer very clear. As we might expect, this transition occurs 
at smaller total mass for higher z. 
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Fig. 3. — Same as Fig. [T] but with mi = 10^ A/0, TO2 = 3 x 10^ Mq, and z = 1 at different values of N (the 
number of days before merger). Reading from left to right, = (solid line), 1 (dashed line), 2 (dash-dotted 
line), 4 (dotted line), 7 (solid line again), 14 (dashed line), 21 (dash-dotted line), and 28 (dotted line). Clearly, 
the largest change — in shape as well as median — happens between merger and one day before. 



The right panel of Figure |3] shows the evolution of the minor axis 26. Again the distribution slowly 
changes shape over time, with the most drastic change occurring in the last day. Here the final distribution 
still retains a slight peak, along with the previously discussed long tail of small errors. Interestingly, this tail 
is present to some degree throughout the evolution. As the total mass increases, the final distribution moves 
to the right until, as with 2a, the sharp change of shape disappears. The same evolution occurs at higher z, 
again with a shift in the mass scale. 

To further understand how the error ellipse evolves, consider Figure |4] Here we show the evolution of 
the median axes for a wide range of masses ("low," "intermediate," and "high") at z = 1. As in Figures [2] 
and [3] our calculations produce output for = 0, 1, 2, 4, 7, 14, 21, and 28; they are connected by lines 
only to guide the eye. Almost all the cases we present show similar behavior: The ellipses gradually shrink 
with time before sharply decreasing in size during the final day. Significant deviations from this behavior 
come from the high-mass binaries, which evolve more drastically at large N before settling in to resemble the 
lower mass curves. This high-mass deviation is an artifact of our choice of maximum N. Smaller binaries 
may spend many months or even years in the LISA band before merger. In those cases, enough signal has 
already been measured at = 28 to locate the binary reasonably well. By contrast, the high-mass binaries 
spend much less time in banc^and have not been measured so well by A'^ = 28. They have to "catch up" 
to the smaller mass binaries over the first few weeks of our measurement window. Nearly identical results 
were found by K07: Figure 2 of K07 plots the evolution of sky position for an intermediate-mass binary from 
A'^ ~ 300. They find that the measurement accuracy rapidly evolves early in the measurement, with slopes 
of angular error versus time very similar to what we show in Figure |4] for high-mass binaries. 

Although these curves are qualitatively quite similar, there are significant quantitative differences. For 
example, the evolution of the intermediate-mass binaries is less than that of the low-mass binaries, especially 



'''The two highest mass binaries in Fig. 4 are not even in band a full 28 days. 
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Fig. 4. — Medians of the sky position ellipse axes for Monte Carlo runs of 10^ binaries as a function of time 
before merger. Major axes 2a are on the left; minor axes 26 are on the right. Data were only outputted at 
the marked points; the lines are there just to guide the eye. The masses have been subdivided into "low," 
"intermediate," and "high" groups; the exact values (in units of solar masses) are given in the legends. Note 
also the different scales (arcminutes and degrees) on the y-axis. 
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in the last day. The intermediate-mass major axes shrink by a factor of 4 — 6 over the entire four-week 
period, whereas the low-mass axes shrink by a factor of ~ 5 — 11. Another important quantitative difference 
can be seen by comparing the major and minor axes. We find that the ratio 2a/2fe grows with time in 
most cases, indicating that the minor axis tends to shrink more rapidly than the major axis. The only 
exceptions are the previously described high-mass cases, in which 2a/2& shrinks during most or all of the 
inspiral. Presumably the same behavior would also be seen for lower mass binaries at higher values of N; 



this conclusion is supported by Figure 4 of Kocsis et al. (2007a). For all other cases, 2a/2b ~ 1.3 — 1.7 at 



= 28 and increases to ~ 1.8 — 4 at = 0. The largest increase is typically in the final day. 

What causes this dramatic improvement in the last day of inspiral? Three factors primarily contribute 
to our ability to localize a source on the sky: modulations due to EISA's orbital motion, modulations due 
to spin precession, and S/N accumulated over time. Since LISA moves the same amount in the final day as 
in any other day, orbital-induced modulations cannot be the cause. A great deal of S/N is accumulated in 
the last day (typically increasing by a factor of ~ 2 or more), and many parameter errors scale as (S/N)~^. 
However, K07 demonstrate that sky position and distance errors do not scale as (S/N)~^ in the last few 
weeks before merger. Our "no precession" code supports their conclusion: The final jump in S/N cannot 
make up for the lack of orbital modulation over such a short timescale. 

The remaining possibility is spin-induced precession. Indeed, we found in Paper I that the number of 
modulations due to spin precession increases dramatically as the binary approaches merger. This suggests 
that the improvement we see is due to the impact of precession. To examine this hypothesis, we plot the 
"precession" and "no precession" results together on the same axes. Figure |5] shows such a plot for a low- 
mass system and an intermediate-mass system at z = 1. We see that at A^ = 28 days, the two codes give 
similar results for localization accuracy. Their predictions gradually diverge as merger is approached. The 
greatest jump between the two codes occurs on the last day before merger, agreeing with our expectation that 
precession effects are maximal then. The effect is greater in the low-mass case than in the intermediate-mass 
case. Similarly, the ratio 2a/26 starts about the same in both codes, growing very slowly until A = 1. At 
this point, it jumps dramatically in the precession code, while staying roughly the same in the no precession 
code. Interestingly, one effect of precession is that the localization errors track (S/N)~^ rather closely. By 
breaking the various correlations which made them deviate from (S/N)^^, precession-induced modulations 
allow the errors to evolve in a manner that is more consistent with our naive expectations. 

The time evolution of the correlations of interest, those between the sky position and binary orientation, 
are illustrated in Figure [g] Here we show the off-diagonal components of the covariance matrix S"'' (where 
a £ {fi-N, 4>n} and b S 1/1^(0), (^l(O)}) for the binary in Figure[2]3. We found that examining the normalized 
correlation coefficients = 5]°''(S'"'S'''')"^/^ can mislead since and E'''' are rapidly evolving at the 
same time as E"''. At large A^, the two sets of angles are relatively strongly correlated due to degeneracies 
in the measured waveform (|6|. However, in the last day before merger, the correlations sharply decrease 
as precession effects accumulate. The reduction of these correlations coincides with the sudden drop in 
parameter errors seen in Figure [2] (and, by extension to the entire Monte Carlo run. Figs. [3]-|5]). 

Finally, we investigate the time evolution of errors in the luminosity distance D]^ . Figure [7] shows the 
distribution of ADl/Dl (determined solely by taking into account GW measurement effects) evolving in 
time for a binary with mi = lO^Af0, TO2 = 3 x 1O^M0, and z = 1. We see that in contrast to sky position, 
the shape of the distribution does not change very much with time. It typically spreads enough to reduce its 
height, but it maintains a well-defined peak. However, the progression of the median is very similar to the 
sky position case: It decreases slowly with time until the last day, when it jumps drastically. The evolution 
of median values of ADl/D^ follows tracks very similar in shape to those shown in Figure [4] so we do not 
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Fig. 5. — Medians of 2a and 2b as a function of time, comparing an analysis that accounts for spin-induced 
precession to one that neglects it. Solid lines trace the evolution of 2a, dashed lines trace 2b. Precession results 
are marked with crosses, no precession with circles. The left plot shows a low-mass case, rrii = 3 x IO^Mq 
and m2 = lO^M©; the right plot shows an intermediate-mass case, mi = 7712 = IO^AIq. Both plots are for 
z=l. 
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Fig. 6. — Off-diagonal covariance matrix entries illustrating correlation between sky position and binary 
orientation as a function of time, for the binary in Fig. ^jp,. The correlations decrease rapidly in the final day 
before merger, when precession effects are maximal. 



show them explicitly. 

Table |5] summarizes all of our results on the time evolution of localization. At low redshift, the ability 
to locate an event on the sky is quite good over much of the mass range even as much as a month in advance 
of the final merger. In most cases, the localization ellipse at z = 1 is never larger than about 10 deg^ in size, 
which is comparable to the field of view of proposed future surveys, such as LSST ( Tyson et al.||2002 ). This 
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Fig. 7.— Distribution of I^Dl/Dl for 10"^ binaries with mi = 10*^ M©, to2 = 3 x IO^Mq, and z = 1 at 
different values of N (the number of days before merger). Reading from left to right, = {solid line), 1 
(dashed line), 2 (dash-dotted line), 4 (dotted line), 7 (solid line again), 14 (dashed line), 21 (dash-dotted line), 
and 28 (dotted line). 



ability degrades fairly rapidly as redshift increases, especially for larger masses. At z 3, the ellipse can be 
~ 10 deg^ a few days in advance of merger for small and intermediate masses. At the highest masses, GWs 
provide very little localization information. Going to z = 5 makes this even worse; an ellipse of ^ 10 deg^ 
can be found at most a day prior to merger, and only for relatively small mass ranges. 

In many cases, ADl/Dl is determined so well by GWs that gravitational lensing errors are expected 
to dominate. As such, the GW-determined values of ADj^/D^ are essentially irrelevant for locating these 
binaries in redshift space; lensing will instead determine how well redshifts can be measured. However, in 
some cases, the intrinsic distance error exceeds the lensing error, so it is worth knowing when the "lensing 
limit" can be achieved. At z = 1, the limit is achieved for most binaries as long as a month before merger. 
For z = 3, the limit is only achieved a few days to a week (depending on mass) before merger; for z = 5, 
intrinsic errors generally exceed the lensing errors even at a day before merger. 



3.3. Results II: Angular dependence of localization accuracy 

We now examine one final interesting property of the errors: their dependence on the sky position of 
the source. As we design future surveys to find counterparts to MBHB coalescences, it will be important 
to understand if there is a bias for good (or bad) localization in certain regions of the sky. It is also useful 
to know in advance whether the "best" regions are likely to be blocked by foreground features such as the 
Galactic center. 

Before discussing this dependence in detail, it is worth reviewing some details of how our Monte Carlo 
distributions are constructed. As described in § |2.4| in most of our analysis we randomly distribute the sky 
position of our binaries, drawing from a uniform distribution in fij^ — cos0n and where On and 
are the polar and azimuthal angles of a binary in solar system barycenter coordinates. We also randomly 
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Table 5: Summary of time evolution of localization accuracy 



Binary Range N (days) 2a (deg) 26 (dcg) AO at (dcg^) /^Dl/Dl 

1 Low (M < 106 Mq) (17-27)^ (6.6 - 13)=^ 0.02-0.07 (2.4 - 4) x 10"^ 

1 (71 - 140)^ (26 - 91)'^ 0.3 - 2.8 (8 - 21) x IQ-^ 

7 (100 - 180)^ (54- 120)^ 0.9-5.1 (1.3-3)xl0-2 

28 (140 -240)^^ (94- 180)^ 2.6-9.6 (2 - 4.3) x 10"^ 

Int. (10^ Mq < M < 4 X 10^ Mq) (31 - 40)^ (9.2 - 22)^^ 0.04-0.18 (3.8-5.6) x lO'^ 

1 (56 -75)*^ (20 -48)*^ 0.17-0.77 (6.4 - 11) x IQ-^ 
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Note. — Typical ranges of sky position and distance measurement accuracy as a function of time until merger for low-, 
intermediate-, and high-mass binaries. Angles are in degrees unless otherwise noted; solid angles are in square degrees. 
"These angles are in arcminutes; all others are in degrees. 

''Some very massive systems are excluded from these data. In those cases, the position and distance are very poorly constrained 
this far in advance of merger. In some cases, the binary is even out of band. 

■^AU of the binaries of this mass and redshift are either very poorly measured or completely out of band this far in advance of 
merger. 
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Fig. 8. — Dependence of the localization errors on /Iat. The major axis 2a of the sky position error ellipse is 
on the left, the minor axis 2b in the center, and the luminosity distance errors l^D^/ on the right. Each 
datum represents the median of 10* binaries with mi = IO^Mq, m2 = 3 x IO^A/q, and z = 1; all other 
parameters are selected randomly (except for Jiiq, whose range is limited to the bin width). 

choose our binaries' final merger time. In this section, rather than distributing the sky position, we examine 
parameter accuracies for particular given positions; all other Monte Carlo parameters are distributed as usual. 
Because we continue to randomly distribute the final merger time, the relative azimuth between a binary's 
sky position and LISA^s orbital position at merger, 6(j) = 4>n — (t>LisA{tc), remains randomly distributed. As 
such, we expect our analysis to efl^ectively average over (jj^, washing out any strong dependence on this angle 
in our analysis. 

We begin by examining the dependence of errors on p,pj. We evenly divide the range — 1 < /Iat < 1 into 
40 bins and run a Monte Carlo simulation with 10* points in each. That is, we pick /Ijv only from the bin 
range, but we pick all other random parameters in the usual manner. The results for a representative binary 
(mi = 1O^M0, 1712 = 3 X lO^Af0, and z — 1) are shown in Figure [s] 

Note that all of the error distributions are symmetrically peaked around the plane of EISA's orbit 
(^Af = 0). Any slight asymmetry is due only to statistical effects. This is reassuring; LISA should not favor 
one hemisphere over the other. There is also additional structure that is parameter dependent. At its peak, 
the major axis 2a is almost 35% greater than the position-averaged median valuc[^of 31.1'. It then decreases 
with IP-nI and reaches a minimum of about 25' for 0.75 < \P-n\ < 0.8. Finally, there are subpeaks near the 
ecliptic poles, although they still lie below the position-averaged median. The dependence of the minor axis 
2b on angle is even more dramatic. At its peak, 2b differs from the position-averaged median of 13.0' by 
over 85%. Just like the major axis, it drops to a minimum, but it does so more rapidly. The minimum also 
occurs at a slightly larger value of \fiN\ than for 2a. The minor axis also shows fairly strong subpeaks near 
the ecliptic poles, with values higher than the position-averaged median. 

The luminosity distance errors behave slightly differently. First of all, the variation with /Iat is weaker 
than for the sky position: At the central peak, ADl/Dl is only ^ 15% greater than its position-averaged 
median value (0.00392). In addition, while the distribution again peaks near the poles, it does so at a smaller 
value of I /I AT I . (In fact, the peaks occur very close to where sky position errors are minimized.) 

By binning the data sets developed for Paper I, we were able to confirm this behavior over a wide range 
of masses and redshifts, albeit with poorer statistics: 10* binaries in total, rather than per bin. Thanks 



"Note that the position-averaged medians quoted here are slightly different from those quoted in Table 1, since this sample 
has 40 times more points. 



- 25 - 






Fig. 9. — Same as Fig. [s] but for dependence on and with 10^ binaries in each bin rather than 10^ 



to the poorer statistics, we cannot resolve the small polar subpeaks in the distribution of 2a. In addition, 
in some cases it appears that the side peaks can be larger than the central peak in the distribution of 
/S.Dl/Dl. Aside from these minor variations, the shapes and relative amplitudes seen in these distributions 
hold robustly over all masses and redshifts we considered. 

We next investigate the dependence of the errors on the azimuthal angle ipiq. To improve the statistics, 
we now calculate 10^ binaries in each bin; the results are shown in Figure |9] The errors have a very weak 
(although nonzero) dependence on ^^r: The maximum deviation from the overall median is only ^ 1% — 2%. 
This is to be expected; as discussed at the beginning of this section, the randomness of our binaries' merger 
times effectively averages over azimuth. If we did not average over azimuth in this way, we would expect 
a moderately strong 0jv-dependence due to the functional form of EISA's response to GWs. Even after 
averaging this dependence away, we might expect some residual (f>N structure due to the "rolling" motion 
of the LISA constellation. The phase associated with LISA's roll angle puts an additional oscillation on the 
measured waves (see the a-dependence in eqs. [47] and [48] of K07, where a encodes the roll angle), and 
we do not average over this angle. Indeed, on close inspection, we can make out roughly two peaks in each 
plot in Figure [9] (although the statistics are still too poor to resolve them clearly) , consistent with the cos 2a 
and sin 2a behavior shown in K07. For comparison, we also examined the (/)Ar-dependence for this mass and 
redshift with precession "turned off." The oscillatory behavior appears very clearly in this case. Because of 
the weakness of the (f)^ behavior, we were unable to easily check it for other masses and redshifts. 

Finally, to give an overall sense as to how localization varies on the sky, we present in Figure [TOjsky maps 
of the median major axis 2a, minor axis 2b, localization ellipse area Afi^r, and distance accuracy ADl/Dl. 
We show these data both in ecliptic and Galactic coordinates. Note that most of the region of small 
error lies outside of the Galactic plane. This potentially bodes well for searches for MBHB electromagnetic 
counterparts — the regions where instruments like LISA "see" most sharply are less likely to be hidden by 
foreground features. Certain portions of the sky that LISA sees well will be easier to search telescopically 
than others. It will be an important task for future surveys over all electromagnetic bands to identify regions 
that are particularly amenable to finding counterparts to MBHB events. 



4. Summary and conclusions 

As discussed at length in Paper I, accounting for the general relativistic precession of the angular mo- 
mentum vectors in an MBHB system has a dramatic impact on what we can learn by observing the system's 
gravitational waves. Spin-induced precession breaks degeneracies among different parameters, making it 
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Fig. 10. — Sky maps of major axis 2a (top row, in arcminutes), minor axis 2b (second row, in arcminutes), 
localization ellipse area Ari^r (third row, in square degrees) and ADl/Dl (bottom row) for LISA observations 
of MBHBs in different parts of the sky. Data in the left column are presented in ecliptic coordinates; data 
on the right are in Galactic coordinates, with the Galactic center at the middle. Note that the level of 
(l)N variation is so small that it would not show up in these figures; accordingly, they are essentially just 
remappings of Fig. [8] 
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possible to measure them more accurately than they could be determined if precession were not present. 
This has a particularly important impact on our ability to locate such a binary on the sky and to determine 
its luminosity distance — the degeneracy between sky angles, distance, and orientation angles is severe in 
the absence of precession. 

Our analysis shows that the improvement that precession imparts to measurement accumulates fairly 
slowly. In using one code which includes the impact of spin precession and a second which neglects this effect, 
we find little difference in the accuracy with which GWs determine sky position and distance for times more 
than a few days in advance of the final merger. The difference between the two codes grows quite rapidly 
in these final days. In the last day alone, the localization ellipse area decreases by a factor of ^ 3 — 10 (up 
to ~ 60 in a few low-mass systems) when precession effects are included. Distance determination is likewise 
improved by factors of ~ 1.5 — 7 in that final day. 

Not all of the precession effects occur in the final days. We saw in Figure [3] that the long tail of small 
minor axes can be seen, to some degree, throughout the inspiral. We could get lucky and find a binary with 
a very small value of 26 weeks before merger. But the improvement in the median that we found in Paper I 
appears to take effect only in the final days of inspiral. Therefore, while precession may in fact help improve 
the final localization of a coalescing binary by a factor of ^ 2 — 10 in each direction, it will not be much help 
in advanced localization of a typical binary. 

Nevertheless, the pixel sizes that we find are small enough that future surveys should not have too 
much trouble searching the region identified by GWs, at least over certain ranges of mass and redshift. 
At z = 1, the GW localization ellipse is ^ 10 deg^ or smaller for most binaries as early as a month in 
advance of merger. (At high masses, the ellipse can be substantially larger than this a month before merger, 
but it shrinks rapidly, reaching a comparable size 1 — 2 weeks before merger.) This bodes well for future 
surveys with large fields of view that are likely to search the GW pixel for counterparts. In addition, GWs 
determine the source luminosity distance so well that the distance errors we find are essentially irrelevant — 
gravitational lensing will dominate the distance error budget for all but the highest masses. 

As redshift increases, the GW pixel rapidly degrades, particularly for the largest masses. Let us adopt 
10 deg^ (the approximate LSST field of view) as a benchmark localization for which counterpart searches 
may be contemplated. At z = 3, this benchmark is reached at merger for almost the entire range of masses 
we considered. As little as a day in advance of merger, however, some of the least massive and most massive 
systems are out of this regime. One week prior to merger, the most massive systems are barely located at 
all (ellipses hundreds of square degrees or larger). The intermediate masses do best, but even in their cases 
the positions are determined with ^ 10 deg^ accuracy no earlier than a few days in advance of merger. The 
resolution degrades further at higher redshift. At 2 = 5, systems with M > 6 x 10^ Mq are not located 
more accurately than ~ 30 deg^ even at merger. Smaller systems are located within ~ 10 deg^ at merger, 
but very few are at this accuracy even one day in advance of merger. The luminosity distance errors also 
increase, so much that they exceed lensing errors a few days to a week before merger at z = 3, and only a 
day before merger at z = 5. This degradation hurts the ability to search for counterparts by redshift and 
subsequently use them as standard candles. 

Our main conclusion is that future surveys are likely to have good advanced knowledge (a few days to 
one month) of the location of MBHB coalescences at low redshift (z ^ 1 — 3), but only a day's notice at 
most at higher redshift (z ~ 5). This conclusion may be excessively pessimistic. As mentioned earlier, recent 
work examining the importance of subleading harmonics of MBHB GWs is finding that including harmonics 
beyond the leading quadrupole has an important effect on the final accuracy of position determination ( |Arun 
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et al. 2007 Trias & Sintes 20081. For most masses, these analyses show a factor of a few improvement in 
position, comparable to the improvement that we find when spin precession is added to the waveform model. 
For high-mass systems, the higher harmonics increase the (previously small) overlap with the LISA band; 
consequently, the improvement can be much larger, up to 2 or 3 orders of magnitude in area. Since these 
two improvements arise from very different physical effects, it is likely that their separate improvements can 
be combined for an overall improvement significantly better than each effect on its own. We plan to test this 
in future work (which is just now getting underway) . 

Finally, we have also studied the sky position dependence of EISA's ability to localize sources. We have 
found that the regions of best localization lie fairly far out of the Galactic plane. However, as emphasized 
by N. Cornish (2007, private communication), a proper anisotropic confusion background might impact this 
dependence. In our calculations, we have assumed an isotropic background, neglecting the likely spatial dis- 
tribution of Galactic binaries. Properly accounting for this background is likely to strengthen our conclusion 
that LISA'S ability to "see" is best for MBHB sources out of the Galactic plane. 
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